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ABSTRACT 

In  the  target  tracking  problem  the  main  purpose  is  to  get  an  accurate  estimation 
of  target  states.  In  order  to  achieve  good  estimation,  the  study  of  target  maneuver 
detection  is  important. 

An  estimator  that  computes  a  constant  input  acceleration  vector  is  first  derived 
in  this  thesis.  It  employs  the  concept  of  least  square  estimation  with  data  consisting 
of  the  residuals  of  the  simple  Kalman  filter.  A  detector  for  sensing  the  target 
maneuver  using  the  input  estimates  based  on  a  fixed  number  of  measurements  is 
developed.  Finally,  the  combination  of  the  estimator,  detector,  and  the  simple 
Kalman  filter  is  developed  to  form  a  tracker  for  maneuvering  targets. 

The  maneuvering  target  tracker  developed  experiences  problems  due  to  errors 
in  maneuver  start-time  detection  and  the  computation  load  associated  with  using  a 
variable-length  window  to  estimate  the  input.  Therefore,  a  modified  input 
estimation  method  for  tracking  a  maneuvering  target  is  presented.  It  uses  only  a 
fixed  number  of  measurements  to  compute  the  input  estimates  and  employs  a 
scaling  factor  to  correct  the  input  estimates  for  feeding  back  to  a  second  Kalman 
filter  (used  for  maneuvering  target  tracking).  The  results  of  target  tracking  using 
these  different  methods  are  presented  in  this  thesis  and  the  modified  input 
estimation  method  turns  out  to  yield  better  tracking  of  maneuvering  targets  in 
some  applications. 
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I.    INTRODUCTION 

Most  tactical  weapons  systems  require  accurate  tracking  of  manned 
maneuverable  vehicles  such  as  aircraft,  ships,  and  submarines.  The  proliferation 
and  increasing  sophistication  of  surveillance  systems  has  generated  a  great  deal  of 
interest  in  algorithms  capable  of  tracking  a  large  numbers  of  objects  using 
measurement  data  from  many  and  possibly  diverse  sensors.  Major  advances  in 
hardware  and  algorithms  have  increased  signal  processing  capabilities  by  one  or 
more  orders  of  magnitude  in  recent  years.  This  has  made  the  measurement  data 
available  for  tracking  even  more  numerous  and  complex,  creating  a  demand  for 
corresponding  advances  in  information  processing  techniques  to  deal  with  them. 

Tracking  is  the  processing  of  measurements  obtained  from  a  target  in  order  to 
generate  an  estimate  of  the  target's  current  state  [  Ref.  1.  ].  In  general,  tracking 
algorithms  require  mathematical  models  of  targets  and  inputs.  There  are  two 
fundamentally  different  types  of  inputs,  random  inputs  and  deterministic  inputs. 
Maneuvering  targets  can  be  characterized  by  an  unknown  input  which  can  be 
modeled  as  a  random  process  or  assumed  to  be  a  nonrandom  disturbance,  or  the 
sum  of  both.  These  unpredictable  changes  in  target  motion  represent  a  major 
challenge  in  tracking  system  design. 

The  usual  technique  applied  to  the  problem  of  tracking  maneuvering  targets  is 
the  Kalman  filter.  In  order  to  achieve  good  tracking  accuracy  for  a  maneuvering 
target,  the  target  input  acceleration  must  be  well  modeled  when  applying  the 
Kalman  filter.  In  general,  the  Kalman  filter  can  yield  good  tracking  for  a 
maneuvering  target  with  an  unknown  random  input  or  a  known  deterministic 
input.  Therefore,  the  application  of  Kalman  filtering  techniques  to  the  problem  of 


tracking  maneuvering  targets  generally  requires  a  stochastic  forcing  function  or  a 
deterministic  input  acceleration.  If  the  forcing  function  is  only  assumed  to  be  a  " 
white  "  Gaussian  process,  the  filter  model  of  target  motion  may  be  unrealistic.  A 
more  realistic  approach  is  to  model  target  motion  as  the  combination  of  random 
input  accelerations  and  deterministic  input  accelerations  or  just  the  deterministic 
input  accelerations.  A  number  of  different  approaches  have  been  applied  to  the 
maneuvering  target  problem  [  Ref.  2.  ]: 

•  the  target  acceleration  was  modeled  as  a  random  process  with  known 
exponential  autocorrelation, 

•  maneuver  is  modeled  as  an  increase  in  the  plant  driving  noise  (  white 
Gaussian  sequence ),  or 

•  the  acceleration  is  assumed  to  be  limited  to  a  time-invariant  set  of  discrete 
values  and  switched  values  according  to  a  semi-Markov  process.  The 
acceleration  is  then  estimated  by  hypothesis  testing. 

There  are  also  many  additional  algorithms  developed  for  tracking  a 
maneuvering  targets.  All  of  the  approaches,  as  mentioned  previously,  assume  some 
a  priori  statistical  description  of  the  maneuver,  ranging  from  white  noise  to 
colored  noise  to  a  semi-Markov  process.  The  algorithms  work  well  within  the 
context  of  their  underlying  assumptions.  If  the  assumptions  made  do  not 
correspond  to  the  actual  nature  of  the  maneuvers,  filter  performance  may  be 
degraded. 

In  this  thesis,  an  algorithm  for  tracking  maneuvering  target  is  presented  and  we 
assume  only  that  the  maneuver  is  a  constant  acceleration.  The  motivation  of  this 
algorithm  came  from  Fig.  1.1  which  shows  a  typical  maneuvering  target  being 
tracked  by  two  different  Kalman  filters.  In  Fig.  1.1a,  the  target  is  tracking  by  a 


simple  Kalman  filter  which  is  a  filter  with  no  prior  knowledge  about  the 
deterministic  input.  The  estimated  state  has  some  bias  relative  to  the  actual  state 
after  maneuver  start.  In  Fig.  1.1b  the  estimated  state  and  the  actual  state  of  the 
target  tracking  using  the  Kalman  filter  with  known  deterministic  input  are  almost 
the  same  before  and  after  the  maneuver  start  and  much  better  than  the  simple 
Kalman  filter. 

The  tracking  algorithm  in  this  thesis  includes  an  algorithm  to  estimate  the 
unknown  deterministic  input  using  the  simple  Kalman  filter  and  to  update  the 
simple  Kalman  filter  using  the  input  estimates.  The  difficulty  here  is  that  the  time 
of  occurrence  of  maneuvers  may  be  unknown.  This  problem  can  be  handled  by 
making  the  filter  simultaneously  perform  both  target  tracking  and  maneuver 
detection.  It  seems  reasonable  to  use  the  observation  residuals  available  in  the 
Kalman  filter  with  statistical  hypothesis  testing  to  detect  the  maneuver.  When  a 
maneuver  is  detected,  the  magnitude  of  the  acceleration  is  identified  using  a  least- 
square  algorithm.  The  result  is  used  in  conjunction  with  a  simple  Kalman  filter  to 
estimate  the  state  of  the  target. 

The  aim  of  the  acceleration  input  estimation  is  to  remove  the  filter  bias  caused 
by  the  target  deviating  from  the  assumed  zero-mean  white-noise  random 
acceleration  motion.  As  is  well  known,  bias  removal  is  accompanied  by  an  increase 
in  the  estimation  variance.  A  detector  is  used  that  checks  the  magnitude  of  the 
estimated  inputs  to  determine  the  occurrence  of  a  maneuver.  The  simple  Kalman 
filter  is  used  alone  during  periods  when  no  maneuver  takes  place.  The  method 
mentioned  above  is  called  the  input  estimation  method. 


Figure  1.1a    Maneuvering  Target  Tracking  Using  the  Simple  Kalman 

Filter 


Figure  1.1b    Maneuvering  Target  Tracking  Using  the  Kalman  Filter 

With  Known  Input 


Combining  the  input  estimation  method  and  the  adaptive  Kalman  filtering 
algorithm,  we  can  detect  the  target  maneuvering.  It  is,  however,  difficult  to  exactly 
detect  the  actual  maneuver  start  time.  Also,  it  is  hard  to  get  good  estimates  of  the 
unknown  deterministic  input  using  a  variable  length  window  of  measurements 
from  the  inaccurate  maneuver  start  time  as  required  by  the  theory.  In  addition,  the 
computation  time  when  using  this  variable  length  window  technique  to  estimate 
input  is  prohibitive.  To  overcome  the  difficulties  in  the  detection  and  the  estimation 
problem,  a  modified  input  estimation  approach  is  proposed  in  this  thesis  which  is 
easier  to  implement.  With  this  new  approach,  the  bias  produced  by  the  maneuvers 
can  be  removed  and  better  tracking  is  achieved  for  maneuvering  targets  when 
compared  with  the  original  method. 

A  brief  introduction  about  the  Kalman  filtering  algorithm  and  target  model  is 
given  in  Chapter  II.  Chapter  III  presents  the  algorithms  for  input  estimation,  state 
correction,  and  maneuver  start  time  detection.  Simulation  results  are  presented  in 
Chapter  IV.  The  conclusions  are  in  Chapter  V. 


II.  TRACKING  THEORY 

A.  TARGET  MODEL 

The  target  model  selected  in  this  thesis  is  based  on  the  assumption  that,  without 
maneuvering,  the  aircraft  generally  has  a  zero-mean,  white-noise  random 
acceleration.  If  the  target  was  not  able  to  deviate  from  this  acceleration  model,  i.e., 
could  not  maneuver,  then  the  tracking  problem  could  be  solved  quickly  and  simply 
using  the  Kalman  filtering  algorithm.  However,  the  maneuver  capability  of  the 
aircraft  can  cause  the  estimated  states  to  have  some  biases  relative  to  the  actual 
states. 

For  the  model  of  the  aircraft  motion,  the  state  equations  are  derived  for  the 
actual  continuous  time  target  motion  and  are  then  expressed  in  discrete  time 
according  to  the  standard  discretization  procedure.  The  model  below  is  presented 
for  three-dimensional  coordinates  in  the  direction  of  range,  elevation  angle,  and 
bearing  angle.  We  will  use  a  single  spatial  variable  later  in  this  thesis  to  simplify 
tracking  performance  evaluation.  This  simplification  is  reasonable  since  the 
estimates  in  each  coordinate  axis  are  nearly  independent.  The  targets  under 
consideration  normally  move  in  response  to  a  zero-mean,  white-noise  random 
acceleration  consisting  of  turns,  evasive  maneuvers  and  accelerations  due  to 
atmospheric  turbulence. 

The  aircraft  equations  of  the  motion  are  modelled  as  follows: 

x(r)  =  Fx(r)  +  Gw(0  (2.la) 

where 
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The  continuous-time  equations  for  the  aircraft  motion  in  each  spatial  variable 
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where  wr  we  and  wb  are  the  zero  mean,  white  noise  accelerations  in  range, 
elevation  angle  and  bearing  angle  respectively.  Combining  these  equations,  the 
continuous  time  state  equations  for  the  aircraft  motion  are 
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The  discrete  time  target  equations  of  motions  are 


(2.3) 


x(k  + 1)  =  <$>(T)x(k)  +  V(r,T)w(k) 


(2.4) 


where  T  is  the  sampling  time  and 
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For  the  model  in  equation  (  2.3  )  it  can  be  verified  that 
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so  the  final  discretized  form  is 
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2r 

0 

0 

T 

0 

0 

T/r 

0 

0 

(2.5a) 
(2.5b) 


0 

0 

w 

T 

II 

we 

2r 
0 

.wb_ 

0 

T/r_ 

(2.6) 


(2.7) 


(2.8) 


where  wt  is  white  Gaussian  noise  sequence  satisfying 


£{wA}  =  0 


(2.9a) 


E{Wyj}  =  Qk8jk  = 


< 

0 

0  " 

0 

< 

0 

0 

0 

»  - 

(2.9b) 


If  the  available  observations  are  range,  elevation  angle  and  bearing  angle,  then 
the  measurement  equation  is  the  following: 


r 

10     0    0    0    0 

e 

= 

0    10    0    0    0 

p_ 

k 

0    0     10    0    0 

lAu 


+ 


(2.10) 


where  vr,  ve  and  vb  are  additive  white  noise  which  is  due  to  the  noisy  measurements 
of  radar  or  beacon  interrogators  and  satisfying 

E{v4}  =  0 

"o-2      0       0 


E{vftvJ}  =  R4^  = 


0      oj      0 

0  0        (T.2 


(2.11a) 


(2.11b) 


When  a  maneuver  occurs,  the  target  model  becomes 

x*+i=<M*+rwt+ru4 


Z*=Hxt+v, 


(2.12a) 
(2.12b) 


where 


0  k    <    n 

u*=i  (2.13) 

„    u  k    >    n 

and  U£  denotes  the  deterministic  input  acceleration  generated  by  a  maneuver  which 

starts  at  k=n.  Note  that  maneuvers  resulting  from  a  constant  input  will  be  addressed 
for  simplicity.  More  general  maneuvers  can  be  included  in  future  work. 

B.    KALMAN  FILTER  EQUATIONS 

Equations  (2.4)  and  (2.12a)  with  <I>(T),  and  T(r,T)  given  by  equations  (2.6)- 
(2.7)  and  uk  by  equation  (2.13)  have  the  form  for  which  the  optimal  linear  filter  is 
identically  the  Kalman  filter.  Other  filters  can,  of  course,  be  used  to  estimate  the 
target  state  vector  xk;  however,  the  Kalman  filter  provides  the  best  performance  in 
terms  of  minimizing  the  mean  square  estimation  error  and  generally  is  easily 
implemented. 

The  Kalman  filter  equations  with  maneuver  input  are: 

X*+l/*  ~  ®Xk/k  +1  U*  (2.14) 

P*+v*=^P*/*°r  +  rQrT  (2.15) 

K*+i  =  p*+i/±H  (Hp*+i/*H7  +  Rj 

■t+yt+i  =  [*■  ~  *m+i"  j*i+i/t  (2. 17) 

x*+v*+i  =  x*+yjt  +  ^-i+i^t+i  —  **xt+i/*  j 


(2.18) 


where 


x 


k+i/k+i  -  minimum  mean  square  estimate  of  xt+]  given  observed  data  up  to 
and  including  time  k+1;  i.e.,  the  filtered  estimate; 
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x*+v*  =  minimum  mean  square  estimate  of  xM  given  observed  data  up  to 

but   not   including   time   k+1;  i.e.,  the  one-sample-ahead 

prediction. 
The  matrices  P*+vt+i,  Pt+Vt  are  the  covariance  matrices  of  the  estimation  error 

and  one- sample-ahead  prediction  error,  respectively.  Initially,  the  calculation  of 
P0/0  must  be  performed  before  using  the  Kalman  filtering  algorithm.  The  Kalman 

filter  is  initialized: 

x0  =  £{x0J  =  0  (2.19) 

The  covariance  matrix  of  the  estimation  error  is  then  initialized: 


0/0 


=  £[(x0-x0)(x0-x0)' 
=  E[xQx0T] 


=  E 


\r    e    b    r    e 


i    b]   ►= 

Jo 


Jo 


of     0    0    0    0    0 

0     cr^     0    0    0    0 
0    0     o\     0    0    0 

0    0    0     o]     0    0 
0    0    0    0     of     0 

0    0    0    0    0     of 


(2.20) 


Equations  (2.14)-(2.17)  are  presented  in  flow  graph  form  in  Figure  2.1  [  Ref. 


1.]. 


The  innovations  sequence  is  defined 


Z*+i  -  Zi+1     Hxt+ 


i/k 


(2.21) 


1  1 


State  at   t^ 
xk 


Transition  to  t 


k+1 


xti=^k+ruk+r^ 


Measurement  at 

t  k+1 

Z^+pHx^i  +  v^-H 


State  estimate  at  f£ 
X  klk 


State  prediction 

i*+//5Eoi*/*+ru* 


Measurement  prediction 
Zk+l/FHxk+l/k 


Innovations 

Z       =Z       — z 

k+1       k+1     k+llk 


Updated  state 
estimate 

x  k+l/k+rx  k+l/lt^lc   Zk 


State  error  covariance  at  t^ 
*klk 


State  prediction  covarianct 

pk+m*pk,fT+*lrT 


Innovation  variance 

S,   =5HP.     7/Hr+R 

k+1       k+llk 


Filter  gain 

Kk+i=i>k+inTs 


-i 

k+A 


Updated  state 
covariance 

pk+ilk+7«-Kk+Ppk+i 


Figure  2.1     Kalman  Filtering  Algorithm 
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and  has  the  important  property  that  it  is  a  zero-mean,  white-noise  sequence  and  is 
Gaussian  if  the  initial  state  x0  and  the  other  noise  processes  are  also  Gaussian.  Its 

covariance  is 


E[zk+iZMT]  =  E[(Zt+1  -Hx4+V4)(Z4+1  -Hx4+1/t)7 


=  HPt+1/tHr  +  R 

=  Sm  (2.22) 

From  Figure  2.1  we  see  that  the  calculation  of  the  error  covariance  matrices  is 
independent  of  the  maneuver  occurrence.  The  innovation  variance  from  the  simple 
Kalman  filter  and  the  one  from  the  Kalman  filter  with  known  input  will  be  the 
same.  In  Chapter  III  we  are  going  to  use  the  innovations  and  innovation  variance  in 
the  algorithm  to  estimate  the  deterministic  input  acceleration. 
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III.    INPUT  ESTIMATION,  DETECTION  AND  STATE 

CORRECTION 

The  basic  problem  posed  by  a  maneuvering  target  is  that  there  exists  a 
mismatch  between  the  modeled  target  dynamics  and  the  actual  target  dynamics. 
Provided  the  target  model  is  correct,  the  Kalman  filter  will  generate  optimum 
estimates  of  the  target  motion.  However,  if  the  target  initiates  and  sustains  a  sudden 
pilot-induced  maneuver,  then  the  assumed  target  model  is  not  correct.  Unless  this 
step  discontinuity  is  accounted  for,  the  Kalman  filter  will  accumulate  errors  and 
possibly  lose  track.  Therefore,  something  more  is  required  in  the  tracking 
equations  to  account  for  the  maneuver  event. 

The  solution  of  this  problem  can  be  divided  up  into  three  different  phases. 
First,  the  maneuver  must  be  detected.  Second,  the  maneuver  must  be  estimated 
based  on  the  simple  Kalman  filter.  Third,  the  simple  Kalman  filter  state  estimate  is 
corrected  to  compensate  for  the  previous  target  maneuver.  These  phases  will  be 
introduced  in  the  following  sections.  Input  estimation  will  be  discussed  first 
because  the  equations  presented  will  be  used  in  detection. 

A.   INPUT  ESTIMATION 

The  Kalman  estimate,  xt+1/Jk+1,  (  from  equations  (2.14)-(2.18)  )  requires 
knowledge  of  uk ,  which  is  the  target  acceleration  input  and  is  not  available  to  the 

filter.  The  target  acceleration  can  be  estimated  and  included  in  the  Kalman  filter  as 
shown  below.  Let  uk  =  0  in  (2.14),  call  this  the  simple  Kalman  filter,  and  denote  its 
estimate  by  xt+1/4+1.  For  simplicity,  let 
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x*+i/*+i  _  x*+i  (3.1a) 


x*+v*+i  —  x*+i  (3  ib) 


Prior  to  time  nT no  maneuvers  occur  so  that  \k=\k.  At  t  =  nT  the  target  starts 
to  maneuver  with  a  sequence  of  inputs  u„,u,,+1, uB+;,where  u  =  u„,u„+p u„+/  is 

a  constant  input  sequence.  The  simple  Kalman  filter  will  continue  to  give  estimates 
according  to  equation  (2.18)  with  uk  =  0 

xk/k  ~  xk/k-\  +  *^k[Zk  -  Hxt/t_,  J 
where 


and 


=  Oxk+Kt+1[Zt+1-HOxt] 

=  (l-Kt+1H)Oxt+Kt+1Zt+] 
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(3.2a) 


xk/k-\  ~  ®*k-yk-i  (3.2b) 

According  to  equation  (3.1a),  we  have 

=  (i-Kku)mk_l  +  Kkzk  (33) 

Defining 

A,  =  (l-KtH)0  (34) 

then  equation  (3.3)  becomes 

x*  =  A^j  +  KkZk 


(3.5) 


=  At+I(AJkxJk_l  +  KtZ4)+Kt+1Zt+1 

=  A*+i A*x*-i  +  A4+1KAZA  4-  Kt+1ZA+1  ^.6) 

and 

x*+2  =  {^-k^^k+i  A*  jxt_i  +  (At+2At+1  j^i^A +  (AA+2  jKt+1ZA+1 

+K/k+2^t+2  (3.7) 

Similarly,  if  the  sequence  u*,uM, u4+;  were  known,  equation  (2.18)  would 

give  the  following  estimates 

xk/k  =xk/k-\  +  KA[Zt  "Hx^.jJ 

where 


(3.8a) 


x*/*-i  —  <Px*-i/*-i  +  1  u*-i  (3.8b) 

So,  using  equation  (3.1b),  we  have 

S*  =  <*>x*-i  +  ruft.j  +  K4[Z4  -  HfOx,.,  +  niw )] 
=  (l-KiH)(Oxt_1+rut_1)+K4Z, 

=  (I  -  KkU)Oxk_,  +  (I  -  K»H)niM  +  Ktz,  (3  9) 

Define 

N,=(i-K,H)r  (310) 

Using  equation  (3.4),  then  equation  (3.9)  becomes 

x*  =  Atx4_!  +  Ntuik_1  +  KtZt  (3.11) 

and 
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'A+l 


=  (I  -  Kt+1H)Ox4  +  (I  -  KtH)rut  +  Ki+1ZA 

=  AA+1  (A.x,.,  +  N»uH  +  KtZt )  +  Nt+1ut  +  KA+1Z, 

=  (At+1At)xt.,  +(AwN>M+(Nw)ut  +  (Aw)KlZl 


'4+1 


+Kt+1Zt+1 


(3.12) 


and 


xt+2  _  (At+2AA+1AAjxA_,  +(At+2At+1NJkJuJt_,  +(AA+2Nt+]  jut 

Therefore,  we  have  two  different  estimates  of  x  (x  and  x)  after  the  occurrence 
of  the  maneuver,  i.e., 

xk  =  {i-Kku)mk_,+Kkzk  (314a) 

xk  =  (I-  K.HjOx,.,  +  (I -  K.Hjru^  +  K,Z,  (3  Mb) 


Now,  define  the  difference  of  those  two  estimates  as  the  following 


A\k  =  \k-\k 


and  recall  that 

x*  =  x* 

up  to  the  instant  the  maneuver  occurred,  k=n. 
This  will  yield 


(3.15) 


(3.16) 


Ax^fl-K^Otx^-x^J  +  fl-K^ru^  (31?) 
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Specifically: 

1.  For  k  <  n,  i.e.,  no  maneuvering,  we  have   xk_x  =  xt_,  and  u^  =0,  so 

Ax,  =  0. 

2.  For  k  =  n,  i.e.,  when  the  maneuver  just  occurred,  we  have  x^  =  1^  ( 
because  n-1  <  n)  and  un-1  =  0  which  is  prior  to  the  maneuver  start,  so  Axt  =  0  is 

still  zero. 

3.  For  k  >  n,  i.e.,  the  maneuver  has  occurred,  we  have  xw-xw  =  Ax^ 
and  uA_j  =  u ,  so 


Ax,  =  (I  -  K,H)<D(x,_1  -  xH )  +  (I  -  K.Hjru^ 
=  (I-KtH)OAxt.1  +  (l-KtH)ru 


=  (l-KtH)(OAxi_1+ru) 


(3.18) 


In  summary,  we  have 
r 


0 


k     <     n 


v. 


Ax,=i  (3.19) 

(i-K.HXOAx^+Tu)  k     >     n 

Next,  we  are  going  to  find  a  general  relation  between  the  estimated  state  with 
the  maneuver  and  the  one  without  the  maneuver.  At  time  tk  =rn+1,  equation  (3.19) 

yields 

Axn+1=(l-Kfl+1H)(OAxfl  +  ru)  (32Q) 

From  equation  (3.19)  we  know  that  Axn  =  0  ,  so  equation  (3.20)  becomes 


Axn+1=(l-Kn+IH)ru 


(3.21) 


Define 


1  8 


then, 


MB+1=I-K„+1H  (322) 


Ax„+1=Mn+1ru  (323) 


From  equations  (3.19)  and  (3.23)  for  k  >  n,  we  have 

Axt  =  (I - K.HXOAx^  +  Hiw )  k>n  (3  24) 

Define  the  general  equation 

Axk=MkTu  (325) 

then,  equation  (3.24)  becomes 

Miru  =  (l-KtH)(OAxt.1+ru) 

MkTu  =  (I  -  K4H)0(Mwru)  +  (I  -  K,H)ru 
hence, 

M,=(l-K1H)(OM,_1+l)  >k>n  (326) 

It  is  obviously  from  equations  (3.19)  and  (3.23)  that  when  there  is  no 
maneuver,  i.e.,  the  time  prior  to  maneuver, 

M*=0  ,k<n  (3.27) 

Therefore,  from  equations  (3.15),  (3.23),  (3.26)  and  (3.27)  we  have 

xk  =  xk  +  MkTu  (328a) 

where 
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0 


k     <     n 


M,=< 


(3.28b) 

(l-K,H)(OMt_1  +  l)  k     >     n 

The  observed  value  of  the  residual  sequence  Z  (with  unknown  maneuver)  can 
be  related  to  the  innovations  residual  sequence  Z  (i.e.,  if  u  were  known)  using 
equation  (3.28) 


zk=zk-mk/k_,  =  zk-nmkX 

Zk=Zk- Hxt/W  =Zk- HOx,.,  - hix,, 

For  different  values  of  k,  we  have 

1 .  For  k  <  n,  because  ut_,  =0,  Zt  =  Zt 

2.  For  k  =  n,  because  uk_x  =  0,  ZA  =  Zt 

3.  For  k>  n>  because  ut_,  =  u , 


(3.29a) 


(3.29b) 


hence, 


Defining 


we  can  have 


z*  -  zk  =  (zk  -  Hmk_,  )-{zk-  Ho^.,  -  hhiw ) 

=  HO(xt_1-xt.1)  +  Hru 

=  H(OMJk_1+l)ru 


Zt=Z4+H(<DMt.I+l)ru 


Bt  =  H(OMt+l)r 


,  k>  n 


(3.30) 


(3.31) 
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zt 


k    <    n 


zt=< 


(3.32) 


I      Zt  +  Bt_,u  k     >     n 

Since  we  are  only  interested  in  estimating  the  input  acceleration  u,  for 
lc  =  n  +  l,n  +  2, ,n  +  s  we  have 


X+1" 

"z„+1" 

"   B„    " 

Z/.+2 

= 

K2 

+ 

Bn+1 

.^n+s. 

Zn+S_ 

_B»+*-l. 

Defining 


X+1" 

A.+2 

_  ^«+*  _ 

£- 


"zn+1" 

Zn+2 

_     *+*_ 

B  = 


then  equation  (3.33)  becomes 


B„    ' 

Hr 

Bn+! 

= 

H(OMn+1+l)r 

Jn«-1. 

H^M^+iJr 

(3.33) 


(3.34) 


Y  =  Bu  +  £  (3.35) 

The  covariance  of  e,  using  the  innovation  covariance  from  the  simple  Kalman 
filter  in  equation  (2.22),  is  given  by 
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z  = 


Sjt+i 

0 

... 

0 

0 

S*+2 

0 

I 

1 

... 

"•. 

0 

0 


0    s 


kis 


(3.36) 


The  problem  now  is  to  estimate  the  unknown  input  acceleration  u  in  equation 
(3.35)  using  the  concept  of  least-square  estimation.  Since  the  elements  of  the 
vector  e  are  independent  and  zero-mean  with  covariance  X,  we  can  find  the 
estimate,  u ,  by  minimizing  the  quadratic  error 


J  =  [Y-Bu]l(I)-,[Y-Bu]  (3.37) 

The  input  acceleration  estimate  that  minimizes  equation  (3.37)  is  obtained  by 
setting  the  gradient  with  respect  to  u  equal  to  zero 


VuJ  =  0 


(3.38) 


yielding 


(3.39) 


2BrIT1[Y-Bu]  =  0 

The  estimated  input  acceleration  is 

U  =  [B'I-'B]-VX-Y  (340) 

Since  the  elements  of  e  are  independent,  zero-mean  random  variables  with 
covariances  £,  the  least-square  estimate  is  unbiased,  i.e.,  substituting  equation 
(3.35)  into  (3.40)  yields 

£[u]  =  [B^I-'b]"1  BrI-I£[Bu  +  e]  =  u 

The  estimation  error  is 


(3.41) 


u  =  u-u 
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=  -[B'X-'B]-Vx-£  (342) 

Thus  the  mean-square  error,  i.e.  the  covariance  of  the  error,  is 

(j2  =  L  =  E[uiiT] 

=  [BTI-1B]"1BrI-1IL-1B[B7'L"1B]"1 

=  [bTi_1b]  (3.43) 

To  calculate  the  input  estimate,  the  previous  s  measurements  must  be  stored  to 
perform  the  calculation  and  the  maneuver  is  assumed  to  occur  at  the  beginning  of 
those  s  measurements.  The  measurements  could  be  used  in  two  ways.  One  is  to 
compute  the  input  estimates  for  maneuver  detection  using  a  fixed-length  window  of 
measurements.  The  other  way  is  to  estimate  the  deterministic  input  from  the 
detected  maneuver  start  time  using  a  variable-length  window  of  measurements.  In 
either  way,  the  information  needed  to  compute  the  input  estimates  are  the  residuals 
from  the  simple  Kalman  filter,  the  covariance  of  these  residuals,  and  the  weighting 
function,  BA,  computed  in  equation  (3.34).  Therefore,  to  estimate  the  input 

acceleration  requires  only  results  from  the  simple  Kalman  filter,  i.e.,  there  is  no 
need  to  have  the  prior  knowledge  of  the  input  acceleration. 

The  input  estimation  algorithm  consists  of  the  estimation  of  the  unknown  input 
(maneuver)  using  the  residuals  from  the  simple  Kalman  filter  starting  when  the 
maneuver  starts.  The  maneuver  onset  time  is  still  unknown.  Methods  of  detecting 
the  maneuver  start  time  will  be  discussed  in  the  following  section. 
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B.    DETECTION 

The  estimation  algorithm  developed  in  the  previous  section  requires 
knowledge  of  the  maneuver  start  time.  Detection  is  the  process  of  detecting  the 
maneuver  start  time.  This  section  develops  the  detection  algorithm. 

The  innovations  and  its  covariance  are 

Z*+i  =  ZA+1  -  HxA+1/t  (3.44a) 

S*+i  =  HPA+]/tH   +  R  (3.44b) 

for  the  nonmaneuvering  case.  On  the  contrary,  when  a  maneuver  occurs  the 
corresponding  innovations  and  its  covariance  are 

z*+i  =  z*+i  ~  H**+i/*  =  z*+i  ~  H^x*+i  ~  Hr  (3.45a) 

Sjt+i  =  HPt+1/tH   +  R  (3.45b) 

Therefore,  as  a  result  of  the  maneuver,  the  residual  sequence  remains  a  white-noise 
process  having  the  same  variance  as  St+1,  but  now  it  has  a  nonzero  mean.  When  a 

maneuver  occurs,  or  if  a  nonlinearity  develops,  it  manifests  itself  as  a  bias  in  the 
residual  sequence.  Thus,  to  detect  the  occurrence  of  the  maneuver,  it  is  necessary 
only  to  monitor  the  bias.  As  long  as  no  bias  develops,  we  continue  to  use  the  zero- 
mean  white-noise  acceleration  tracker;  when  a  bias  is  detected,  then  we  switch  to 
the  input  estimation  tracker. 

Mathematically,  detecting  the  bias  reduces  to  the  following  hypothesis  test: 

HQ  :  No  Maneuver  Occurs  :  Zk  =  Zk 

//,  :  Maneuver  Occurs  \Zk  =  Zk  -B^u  (3.46) 

where  Zk  is  zero-mean,  white  noise  with  variance  Sk. 
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Detecting  a  maneuver  is  equivalent  to  detecting  the  presence  of  a  deterministic 
signal  of  unknown  amplitude  and  time  of  arrival  in  a  background  of  zero-mean, 
white  noise.  Since  the  measurement  noises  are  assumed  to  be  Gaussian,  then  the 
optimum  test  for  deciding  upon  H0  or  Hl  is  the  generalized  likelihood  ratio  test. 

The  estimation  of  the  input  is  done  according  to  equation  (3.35)  using  a  fixed 
number  of  measurements  for  the  detector.  An  estimate  is  accepted,  i.e.,  a  maneuver 
is  declared  detected,  only  if  it  is  "statistically  significant".  Therefore,  we  must  find 
a  test  statistic  and  compare  it  with  a  given  threshold  to  see  if  there  is  a  maneuver  has 
occurred.  In  equations  (3.36)  and  (3.40)  the  estimated  input  is  unbiased  with 
covariance  a2.  Using  the  theorem  of  hypothesis  testing  with  significance  [  Ref.  3.  ], 
consider  the  significance  test  for  the  estimate  u  is 

\u\>TH  (347) 

where  TH  is  a  threshold. 

The  threshold  is  selected  by  noting  that  if  the  input  is  zero  (  u  =  0  ),  then 
the  mean  of  u  is  also  zero  and  u  should  have  the  following  density  function 

u~N(0,cr2)  (348) 

and  TH  is  chosen  such  that  the  probability  of  false  alarm  is 

p{\u\>TH}=a  (349) 

with  a  =  10"2  or  smaller. 

The  value  of  TH  and  the  probability  of  detection  can  be  obtained  from  tables, 
given  the  probability  of  false  alarm,  p  .  The  performance  of  this  detector  is 

characterized  by  the  probability  that  a  maneuver  is  correctly  detected  when  it 
occurs,  PD,  and  the  probability  that  it  erroneously  declares  a  maneuver  when,  in 
fact,  no  maneuver  was  undertaken,  PFA.  For  example,  if  the  probability  of  false 
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alarm  is  0.01%,  i.e.,  confidence  level  is  99.99%,  then  the  threshold  value  can  be 
determined  from  the  table  of  normal  distribution  function  with  such  probability  of 
false  alarm.  The  threshold  value  corresponds  to  PFA  =0.01%  is  3.89a.  Therefore, 
as  soon  as  the  input  estimate  is  calculated,  its  absolute  value  is  compared  with  the 
threshold  to  see  whether  a  maneuver  has  occurred. 

The  calculation  of  the  input  estimate  u  assumes  that  the  maneuver  starts  s 
measurements  ago.  It  is  necessary  to  go  back  and  compute  the  input  estimates  based 
on  a  variable  length  of  measurements  and  correct  the  state  estimate  of  the  simple 
Kalman  filter  using  the  estimated  input  once  the  maneuver  is  detected.  The  detected 
maneuver  start  time  may  not  be  the  same  as  the  actual  maneuver  start  time.  Because 
the  bias  is  monotonically  increasing  with  time,  the  estimated  input  does  not  have  a 
sudden  jump  at  the  time  of  maneuver  start.  The  way  to  overcome  this  problem  is  to 
find  a  relation  between  the  delayed  maneuver  start  time  and  the  estimated  input. 
We  can  use  this  relation  to  yield  a  mapping  to  correct  the  detected  start  time. 
Although  the  detector  is  not  optimal,  it  should  perform  well  in  a  practical 
operational  environment. 

C.  STATE  AND  ERROR  COVARIANCE  CORRECTION 

When  a  maneuver  is  detected,  the  state  and  error  co variance  must  be  corrected 
as  follows  [  Ref.  3.  ]. 

1.     State  Correction 

All  the  information  gained  from  the  previous  sections  is  now  available  to 
correct  for  the  target's  maneuver.  Recall  from  equation  (3.28)  that 

xt=x,  +  MJu 

where 
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0  k     <     n 


Mk=< 


<      (l-KtH)(OM4.1  +  l)  k     >     n 

So,  the  simple  Kalman  filter  estimate  is  corrected  or  updated  by  u 
through  the  equation 

xVJ+,  =  xk+s+l+MkTu  (3>50) 

where  x"*+*+i  denotes  the  updated  estimate  at  time  (k+s+l)T. 

Since  u  is  unbiased,  taking  the  expectation  in  equation  (3.40)  shows  that 

£[x"t+rti]  =  £[xi+s+1]  (3.51) 

Thus,  the  correction  will  remove  the  bias  in  xt+J+1 . 
2.     Error  Covariance  Correction 

The  bias  removal  is  accompanied  by  an  increase  in  the  estimation 
variance.  A  new  covariance  estimate  must  be  computed  to  reflect  the  change  in 
equation  (3.40).  Subtracting  the  correct  value  from  the  estimated  value,  i.e., 
subtracting  equation  (3.40)  from  (3.28),  we  find  that 

xL+i=x*+,+i+M<r(u-u) 
where  the  assumption  uk+l  =  u  for  /  =  l  to  s  is  used. 

Recall  from  equations  (3.35),  (3.40)  and  (3.43)  that 

u  =  [Br2T1B]"1BTE-IY 


(3.52) 


=  LB'I  -'(Bu  +  £) 
^LB^Bu  +  LB7!-^ 
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=  u  +  LB7I  ~le 

Hence, 


u-u  =  LB'Z_1e  (3.53) 

The  input  estimation  error,  u-u,  is  uncorrected  with  xt+J+1-xt+J+, 

because  the  innovations  form  a  white  noise  sequence.  Thus,  the  covariance  of  x£+,+i 

is 

=  P^+1  +  (Mir)L(M,r)r  054) 

where  (M4r)L(Mtr)  represents  the  increase  in  variance  as  a  result  of  updating 
xt+rt.j  using  u.  Whenever  updating  occurs,  (M^rJl^M^r)  must  be  added  to  the 
value  of  P  in  equation  (2.16)  causing  a  corresponding  increase  in  the  Kalman  gain. 

A  maneuver  is  considered  finished  when  the  input  estimate  based  on 
measurements  from  the  sliding  window  of  fixed  length  s  becomes  insignificant. 
The  length  s  is  a  design  parameter,  as  mentioned  in  Section  B.  We  need  to  find  a 
appropriate  value  of  s  in  order  to  produce  a  reliable  estimate.  We  are  going  to 
examine  this  issue  by  using  a  simulated  maneuvering  target  model  in  the  following 
chapter. 
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IV.     IMPLEMENTATION  AND  SIMULATION  RESULTS 

The  integration  of  the  estimator,  detector,  and  the  simple  Kalman  filter  into  a 
tracking  scheme  is  straightforward.  The  simple  Kalman  filter  is  first  put  into 
operation  and  the  estimator  and  detector  are  activated  when  the  filter  is  near 
steady-state.  In  this  chapter,  we  generate  a  maneuvering  target  model  and  use  this 
model  to  exercise  the  input  estimation  and  detection  algorithms  which  were 
mentioned  in  the  previous  chapters. 

A.  THE  TARGET  MODEL 

The  target  model  equations  are  given  in  Chapter  II.  The  parameters  T  =  1  sec, 
(2  =  0.1,  R  =  2  degree 2  are  used  in  these  equations.  The  initial  conditions  of  the 
estimate  are  taken  to  be  the  first  measurement,  i.e., 

x(0)  =  Z(0)  degree,  x(0)  =  0 
The  formula  for  the  initial  covariance  matrix  P(0/0)  is  given  in  Chapter  II, 
with  the  given  variance  of  the  elevation  angle  (  a]  =  2  )  and  the  variance  of 
derivative  of  the  elevation  angle  (  o]  =  3  ) 

P(0/0)  =  [2    0    ;    0    3] 

The  target  starts  a  constant  acceleration  at  T  =  300  seconds  and  maintains  this 
maneuver  till  7  =  600  seconds  Using  the  MATLAB  program  elevsim.m,  we  can 
simulate  this  maneuvering  target.  Figures  4.1a  and  4.1b  show  the  state  estimates 
(the  elevation  angle  and  the  derivative  of  elevation  angle)  from  the  simple  Kalman 
filter  versus  the  states  of  the  actual  target.  We  see  that  the  estimated  states  have  a 
bias  relative  to  the  actual  states  after  the  maneuver  starts.  The  state  estimates 
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Figure  4.1a    Position  Tracking  Using  the  Simple  Kalman  Filter 
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Figure  4.1b    Velocity  Tracking  Using  the  Simple  Kalman  Filter 
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from  the  Kalman  filter  with  known  input  versus  the  actual  target  states  are  shown 
in  Figures  4.2a  and  4.2b.  We  see  that  the  estimated  states,  using  the  Kalman  filter 
with  known  deterministic  input,  and  the  actual  states  are  almost  the  same  before 
and  after  the  maneuver  starts.  Note  that  the  estimated  states  in  Figure  4.2  are  much 
better  than  the  estimated  states  in  Figure  4.1.  Figures  4.3  and  4.4  show  the  error 
covariance  and  the  filter  gain  for  each  time  step  for  both  of  the  Kalman  filters. 
Figures  4.3  and  4.4  demonstrate  that  the  time  for  both  Kalman  filters  to  converge 
to  steady-state  is  the  same  and  both  have  the  same  error  covariance.  From  the 
above  figures,  we  can  easily  see  that  the  estimated  states  of  the  two  different 
Kalman  filters  are  greatly  different  due  to  the  maneuver.  In  the  following  sections, 
we  are  going  to  use  this  example  to  perform  the  input  estimation  and  input 
detection  algorithms. 

B.   ESTIMATOR  AND  DETECTOR 

In  applying  the  input  estimator  and  detector  to  a  tracking  scheme,  two  factors 
need  to  be  considered.  The  number  of  measurements  s  used  to  estimate  u  for  the 
detector  and  the  value  of  the  threshold  (TH)  in  the  detector  [  Ref.  3.  ].  We  will 
explain  this  statement  in  the  following  sections. 

1.     Determine  the  number  of  measurements  to  compute  input 
estimates 

When  computing  the  input  estimate,  it  is  assumed  that  the  input 
acceleration  is  constant  for  a  period  of  time.  Therefore,  the  input  can  be  detected 
using  a  sliding  fixed-length  measurement  window  of  length  s.  To  assess  the 
performance  of  the  least  squares  estimator  and  to  provide  a  guide  for  selecting  s, 
the  following  experiment  was  performed.  A  nonmaneuvering  target,  i.e.,  u  =  0,  is 
tracked     by     a     simple     Kalman    filter    and    the    input    estimate,     u, 
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Figure  4.2b    Velocity  Tracking  Using  the  Kalman  Filter  with  Known 

Input 
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Figure  4.3a    Error  Covariance  Pll  from  the  Simple  Kalman  Filter 
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Figure  4.3b    Error  Covariance  P12  from  the  Simple  Kalman  Filter 
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Figure  4.3c    Error  Covariance  P21  from  the  Simple  Kalman  Filter 
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Figure  4.3d    Error  Covariance  P22  from  the  Simple  Kalman  Filter 
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Figure  4.4a    Filter  Gain  Kll  from  the  Simple  Kalman  Filter 
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Figure  4.4b    Filter  Gain  K21  from  the  Simple  Kalman  Filter 
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is  computed  from  equation  (3.40)  for  different  values  of  s.  Figures  4.5a  and  4.5b 
show  the  estimated  input  for  the  nonmaneuvering  target  for  two  different  values  of 
s.  For  the  larger  value  of  s,  the  variation  of  the  estimated  input  is  reduced  when 
compared  with  the  smaller  value  of  s.  Table  I  summarizes  the  statistics  of  the 
estimated  inputs  for  different  values  of  s.  As  expected,  the  standard  deviation  of 
input  estimates  (  uest )  decreases  with  increasing  values  of  s.  Table  I  can  be  used  to 
select  a  suitable  value  of  s.  Figures  4.6a  and  4.6b  also  illustrate  the  effects  of 
different  values  of  s  on  the  input  estimates  for  the  maneuvering  target.  They  show 
similar  results  to  those  in  Figures  4.5a  and  4.5b.  The  choice  of  s  determines  the 
accuracy  of  the  least  squares  estimate,  i.e.,  the  accuracy  increases  with  the  number 
of  measurements.  However,  large  values  of  s  will  create  the  problem  of  a  long 

memory  in  the  estimator  and  reduce  the  importance  of  the  most  recent 
measurements.  Therefore,  we  use  s  equal  to  50  measurements  for  this  thesis  in  the 
simulation  runs.  The  program  used  to  generate  these  results  is  elevdetecl.m. 
2.     Determine  the  threshold  (  TH  ) 

From  Table  I,  we  see  that  the  estimated  input  u  has  a  standard  deviation 
equal  to  0.9183  when  using  50  measurements.  The  threshold  value  can  be 
determined  from  this  number.  We  choose  the  TH  to  be  greater  than  3.89  times  the 
standard  deviation,  i.e., 

77/ >  4.24x0.9183  =  3.893 

This  choice  of  threshold  was  found  by  setting  the  false  alarm  rate  smaller 

than  0.01%.  The  next  step  is  to  perform  the  detection  algorithm  with  this  threshold 

value.  Figure  4.7  shows  the  result  that  for  this  particular  target  with  50 

measurements  to  compute  input  estimates,  we  can  detect  that  there  is  a 
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Figure  4.5a  U  Estimates  of  the  Nonmaneuvering  Target  (  s=30  ) 
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Figure  4.5b  U  Estimates  of  the  Nonmaneuvering  Target  (  s=50  ) 
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Table  I    Statics  of  Input  Estimates  for  the  Nonmaneuvering  Targets 
Using  Different  Numbers  Of  Measurements 


max(uest) 

mean(uest) 

std(uest) 

s=30 

6.9278 

-0.8962 

2.4031 

s=40 

4.8421 

-0.5148 

1.4231 

s=50 

3.4207 

-0.3095 

0.9183 

s=60 

1.5392 

-0.2070 

0.6228 

s=70 

1.2294 

-0.1581 

0.5254 

s=80 

1.0435 

-0.1338 

0.4336 

s=90 

0.9133 

-0.1214 

0.3629 

s=100 

0.8153 

-0.1119 

0.3132 
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Figure  4.6a    Input  Estimates  of  the  Maneuvering  Target  (  30  samples  ) 
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Figure  4.6b    Input  Estimates  of  the  Maneuvering  Target  (  50  samples  ) 
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Figure  4.7    Maneuver  Detection  at  K=363  seconds 
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Figure  4.8    U  Estimates  Using  the  Variable  Length  Window 
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maneuver  starting  at  k=363   seconds.    This   result   comes    from   program 
elevdetec.m. 

C.   DETERMINE  ERRORS  IN  DETECTED  START  TIME 

Although  we  have  detected  the  maneuver  start  time  at  k=363  seconds  in  the 
example  displayed  in  Figure  4.7,  it  is  not  the  true  maneuver  start  time.  Errors  in 
the  maneuver  start  time  will  cause  errors  in  the  estimated  input.  In  order  to  find 
this  error  we  can  apply  a  variable  length  window  to  compute  the  input  estimates 
once  the  maneuver  starts.  We  assume  a  false  maneuver  start  time  several  time  steps 
later  than  the  actual  maneuver  start  time  and  then  perform  the  simple  Kalman  filter 
and  start  the  input  estimation  algorithm  with  a  variable-length  window.  The  result 
for  a  false  maneuver  start  time  at  k=315  seconds  is  shown  in  Figure  4.8  and  its 
sample  mean  value  (averaged  from  k=399  to  k=599  seconds)  is  equal  to  3.2131. 
The  large  values  of  the  input  estimates  between  k=315  to  k=400  seconds,  in  Figure 
4.8,  are  due  to  the  short  length  window  used  in  the  beginning  of  the  input 
estimation  computation.  Better  estimates  are  obtained  when  the  window  becomes 
long  enough.  Table  II  summarizes  the  results  taken  over  20  computer  runs  for 
different  false  maneuver  start  time.  Those  results  are  from  program 
IE_sensitive.m.  From  this  table  we  can  determine  a  mapping  to  correct  the  u 
estimate  for  errors  in  detected  maneuver  start  times.  This  corrected  input  estimate 
can  be  used  to  update  the  states  of  the  simple  Kalman  filter. 
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Table  II    U  Estimates  Using  the  Variable  Length  Window  For 
Different  Maneuver  False  Start  Times 


Different  False  Maneuver  Star 

t  Time 

Run# 

MPF=300 

MPF=305 

MPF=310 

MPF=315 

MPF=320 

MPF=325 

1 

2.9590 

3.0998 

3.2660 

3.3627 

3.5420 

3.7660 

2 

2.9568 

3.0303 

3.1897 

3.3630 

3.4558 

3.6601 

3 

3.0268 

3.1122 

3.3504 

3.4380 

3.5618 

3.8521 

4 

2.9180 

3.1625 

3.2000 

3.3178 

3.5992 

3.6722 

5 

3.0367 

3.0629 

3.2577 

3.4503 

3.5023 

3.7537 

6 

2.9973 

3.0807 

3.3051 

3.3983 

3.5190 

3.8037 

7 

3.0144 

3.1770 

3.3063 

3.4228 

3.6264 

3.8096 

8 

2.9685 

3.1289 

3.2284 

3.3757 

3.5702 

3.7076 

9 

3.0496 

3.1213 

3.2813 

3.4672 

3.5679 

3.7769 

10 

3.0340 

3.2192 

3.2646 

3.4504 

3.6770 

3.7500 

11 

2.9951 

3.1521 

3.2131 

3.4022 

3.6081 

3.6991 

12 

2.9316 

3.0964 

3.3314 

3.3342 

3.5385 

3.8347 

13 

2.9868 

3.1348 

3.2805 

3.3939 

3.5818 

3.7733 

14 

2.9515 

3.0843 

3.2378 

3.3495 

3.5190 

3.7153 

15 

2.9996 

3.1278 

3.2196 

3.4019 

3.5724 

3.6994 

16 

3.0137 

3.0707 

3.2660 

3.4369 

3.5174 

3.7530 

17 

2.9922 

3.0679 

3.2980 

3.3932 

3.5009 

3.7916 

18 

2.9685 

3.0016 

3.2255 

3.3704 

3.4239 

3.7074 

19 

3.0050 

3.0930 

3.2545 

3.4102 

3.5360 

3.7478 

20 

3.0067 

3.1520 

3.2371 

3.4119 

3.6019 

3.7219 

mean 

2.9911 

3.1088 

3.2607 

3.3975 

3.5511 

3.7498 

std. 

0.0359 

0.0512 

0.0432 

0.0404 

0.0588 

0.0524 
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D.  FINAL  RESULTS 

Once  the  maneuver  start  time  is  detected,  we  are  going  to  update  the  state  and 
error  covariance  as  given  in  Chapter  III,  Section  C.  In  addition  we  have  tried 
several  different  approaches  suggested  in  the  course  of  this  research. 

1.  Method  A 

This  method  is  based  on  the  original  idea  mentioned  in  Chapter  III.  It 
updates  the  state  and  error  covariance  as  given  in  Chapter  in  as  soon  as  the  absolute 
value  of  the  estimated  input  becomes  greater  than  the  threshold  value.  Program 
elevesdcr.m  performed  this  method.  Figure  4.9  shows  the  target  track  using  this 
method  and  Fig.  4.10  shows  the  square  difference  of  the  estimated  target  states 
from  the  Kalman  filter  with  known  input  and  the  estimated  states  from  method  A. 
As  Figures  4.9  and  4.10  show,  method  A  apparently  is  not  effective  in  keeping  the 
position  and  velocity  errors  small.  On  checking,  it  soon  became  apparent  that  this 
was  because  the  residuals  from  the  Kalman  filter  including  input  correction  are 
used  in  the  input  estimation  algorithm.  In  turn,  this  produces  transients  and  large 
errors.  This  led  us  to  revise  our  strategy  slightly. 

2.  Method  B 

This  method  tried  a  different  way  to  update  the  state  and  error 
covariance.  It  updated  the  state  and  error  covariance  based  on  Chapter  III  only 
when  the  absolute  values  of  the  estimated  input  is  not  greater  than  the  threshold 
value.  Program  elevesdcr_a.m  performed  this  algorithm.  Figure  4.11  and  4.12 
show  the  result.  As  we  can  see,  the  errors  produced  by  this  method  are  quite  small 
comparative  to  the  results  of  method  A. 
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3.     Method  C 

Since  it  is  hard  to  detect  the  true  maneuver  start  time,  we  can  use  a  new 

scheme  to  solve  the  tracking  problem.  In  this  method  we  track  the  target  using  a 

simple  Kalman  filter  to  generate  the  estimated  input  using  a  fixed  number  of 

measurements  and  then  feed  this  estimated  input  back  to  a  second  Kalman  filter. 

The  results  are  shown  in  Figures  4.13  and  4.14  which  is  not  a  very  good  tracking. 

The  reason  is  shown  in  Figures  4.15,  we  see  that  the  estimated  input  has  a  transient 

that  lasts  for  about  300  samples  and  the  final  value  is  quite  large  compared  to  the 

true  input  (  u=3g  ).  Therefore,  we  need  a  factor  to  scale  down  the  estimated  input 

and  Figure  4.16  shows  the  result  for  a  factor  equal  to 

magnitude   of   true   input 

factor  =  — ~     ~ ^— f— 

mean(uest(600:999)) 

(in  this  example,  the  factor  is  0.1209).  The  scaled  input  estimates  are  shown  in  Fig. 
4.17.  The  squared  position  and  velocity  error  are  shown  in  Fig.  4.18.  Comparing 
Figures  4.16  to  4.18  with  Figures  4.9  to  4.12,  we  see  that  using  the  scaled  input 
estimates  in  the  Kalman  filter  yields  better  tracking.  This  method  provides  an 
easier  way  to  solve  the  tracking  problem.  We  can  find  a  value  of  the  scaling  factor 
as  a  function  of  the  number  of  measurements.  This  method  has  biases  that  last  for 
about  300  time  steps,  then  the  estimates  will  be  close  to  the  actual  states.  Table  III 
shows  the  scaling  factor  for  different  magnitudes  of  true  input  and  for  different 
numbers  of  measurements.From  Table  III,  it  is  easily  seen  that  the  scaling  factor 
only  depends  on  the  number  of  measurements  which  were  used  to  compute  the 
input  estimates  (not  on  the  input  magnitude).  Therefore,  we  can  use  this  table  to 
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determine  the  scaling  factor  for  the  estimated  input  based  on  the  number  of 
measurements  and  apply  the  scaled  input  estimates  to  the  simple  Kalman  filter  to 
yield  better  tracking  accuracy. 


Table  III    The  Scaling  Factors  For  The  Estimated  Inputs 


Magnitude   of  the   true   input 

Measurements 

u=2g 

u=3g 

u=4g 

s=50 

0.1205 

0.1209 

0.1202 

s=60 

0.1615 

0.1621 

0.1611 

s=70 

0.2049 

0.2056 

0.2043 
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V.  CONCLUSION 

In  this  thesis,  a  method  of  tracking  a  maneuvering  target  has  been  presented. 
This  scheme  incorporates  a  simple  Kalman  filter,  an  input  acceleration  estimator, 
and  a  detector.  In  order  to  get  accurate  tracking  of  the  maneuvering  target,  we 
need  to  detect  the  maneuver  start  time  and  apply  the  estimated  input  to  the  simple 
Kalman  filter.  The  input  estimator  presented  in  this  thesis  generates  good  estimates 
of  the  unknown  input  acceleration  when  using  a  sufficient  number  of 
measurements.  The  detection  of  the  maneuver  start  time  is  difficult;  it  is  always 
biased  by  a  time  delay  and  most  of  the  time  spent  on  this  thesis  was  spent  on 
determining  an  accurate  start  time.  We  can  use  the  variable-window  to  correct  this 
time  delay  bias.  Once  we  find  the  correction  factor  based  on  Table  II,  the  estimated 
input  can  be  updated  to  approach  the  true  input  and  the  maneuver  start  time  could 
be  reset.  The  detection  algorithm  in  this  scheme  requires  a  significant  amount  of 
computation  and  memory,  although  it  is  simple  in  concept.  The  advantage  of  the 
detection  algorithm  is  that  it  is  implementable  without  any  a  prior  knowledge  of 
the  maneuvering  characteristics  of  the  target. 

Since  it  is  hard  to  detect  the  exact  maneuver  start  time,  we  have  modified  this 
tracking  scheme  slightly.  We  use  the  simple  Kalman  filter  to  generate  the  input 
estimates  using  a  data  window  of  fixed  length.  In  order  to  produce  reliable  tracking 
of  the  maneuvering  target,  the  estimated  input  must  be  scaled  to  the  true  magnitude 
of  the  input  acceleration.  Fortunately,  Table  III  provides  us  a  way  to  do  this.  The 
scaling  factor  can  be  easily  determined  and  depends  only  on  the  number  of 
measurements  which  are  used  to  compute  the  input  estimate.  Therefore,  as  shown 
in  Fig.  4.16,  this  new  method  can  yield  good  tracking  accuracy  for  the  constant 
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acceleration  considered  in  this  thesis.  There  still  exists  some  biases  at  the  time  of 
maneuver  initiation,  but  this  dies  out  eventually.  Since  we  use  a  one  second 
sampling  time  in  this  thesis,  the  time  for  the  biases  to  be  vanished  is  about  3  minutes 
which  is  not  practical  in  some  applications.  However,  this  method  can  be  applied  in 
applications  with  higher  sampling  rates. 

The  analysis  in  this  thesis  assumed  a  very  specific  situation,  viz.,  a  two-state 
one-dimensional  elevation  angle  tracker.  A  future  step  is  to  incorporate  the 
methods  mentioned  in  this  thesis  into  a  realistic  three-dimensional  tracker  as 
mentioned  in  Chapter  II.  Also,  it  is  important  to  generalize  these  methods  to 
account  for  different  maneuver  models.  Further  reductions  in  the  computation 
load  are  possible  if  the  coefficients  Mk  can  be  assumed  to  be  time-invariant  or 
otherwise  vary  slowly  with  time,  allowing  the  Mk's  to  be  precalculated.  This  would 
greatly  simplify  the  required  real-time  processing  load  [  Ref.  4.  ]. 
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APPENDIX 


COMPUTER  PROGRAMS 


1.  ELEVSIM.M 

2.  ELEVTRACK.M 

3.  ELEVDETECI.M 

4.  ELEVDETEC.M 

5.  EE_SENSITIVE.M 

6.  ELEVESDCR.M 

7.  ELEVESDCR_A.M 

8.  ELEVESDCR_D.M 

9.  ELEVESDCR  Dl.M 


Maneuvering  and  nonmaneuvering  target,  simulation 

Uses  the  simple  Kalman  filter  and  the  Kalman  filter 

with  known  input  to  track  the  target. 

Computes  the  input  estimates  for  the  different  targets. 

Performs  the  detection  algorithm. 

Computes  the  errors  for  false  maneuver  start  times. 

Tracking  solved  by  Method  A. 

Tracking  solved  by  Method  B. 

Tracking  solved  by  Method  C. 

Sub-program  of  ELEVESDCR_D.M 
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1.    ELEVSIM.M 

%  File  Name  :  elevsim.m 

% 

%  Purposes  : 

%     1.  Simulate  the  unmaneuvering  and  maneuvering  target 

%    2.  Elevation  angle  only. 

% 

%  By  :  Meng,  Hsing-Han 

clear 

clc,disp('         '); 

fname=input('Give  the  name  for  the  target  ?  ','s'); 

N=input('Give  the  total  number  of  data  ?  '); 

MP=input('Give  the  time  to  start  maneuver  : '); 

amp=input('Give  the  magnitude  of  input : '); 

%  TARGET  SIMULATION 

rand('normar); 

n=l:N; 

var_we=0. 1 ; 

var_e=2; 

w=l  *sqrt(var_we)*rand(n); 

v=sqrt(var_e)*rand(n); 

u=[zeros(l,MP-l)amp*ones(l,N-(MP-l))]; 

T=l; 

phi=[lT;01]; 
H=[l  0]; 
x(:,l)=[23;0]; 

z(l)=H*x(:,l)+v(l); 

x_manu(:,l)=[23;0]; 
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z_manu(  1  )=H*x_manu(:,  1  )+v(l ); 
gam=[(TA2)/(2*  1 000);T/1 000] ; 
gam_manu=gam; 

fork=l:N-l 

x(:,k+l  )=phi*x(:,k)+gam*w(k); 
z(k+l)=H*x(:,k+l)+v(k+l); 

x_manu(:,k+l)=phi*x_manu(:,k)+gam_manu*(w(k)+u(k)); 
z_manu(k+l)=H*x_manu(:,k+l)+v(k+l); 
end 

%  PLOTTING 

i=l:N; 
for  j= 1:2 

clg 

subplot(21  l),plot(i,x(j,:)),grid 

subplot(212),plot(i,x_manu(j,:)),grid, pause 
end 
clg 

subplot(21  l),plot(i,z),grid 
subplot(212),plot(i,z_manu),  grid,  pause 
clg 

%  SAVE  TO  WORKSPACE 

eval(['save  ',fname,'_u  x  z  N  MP']) 

eval(['save  ',fname,'_m  x_manu  z_manu  N  MP']) 
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2.    ELEVTRACK.M 

%  File  Name  :  elevtrack.m 

% 

%  Purposes  : 

%     1 .  Track  the  unmaneuvering  and  maneuvering  target 

%    2.  Elevation  angle  only. 

% 

%  By  :  Meng,  Hsing-Han 

clear 

clc,disp('        '); 

fname=input('Give  the  name  of  target  to  be  processed  ?  ','s'); 

amp=input('Give  the  magnitude  of  input  : '); 

eval(['load  ',fname,'_m']); 

x=x_manu; 

z=z_manu; 

%  INITIAL  CONDITION 

var_we=0.1; 

var_e=2; 

var_ed=3; 

u=[zeros(  1  ,MP- 1 )  amp*ones(  1  ,N-(MP- 1 ))] ; 

T=l; 

phi=[lT;01]; 

H=[l  0]; 

gam=[(TA2)/(2*  1 000);T/(  1 000)] ; 

I=eye(2); 

Q=var_we; 

R=var_e; 

xe(:,l)=[z(l);0]; 

xem(:,l)=[z(l);0]; 

P=diag([var_e  var_ed]); 
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Pm=diag([var_e  var_ed]); 
i=0; 

%  TARGET  TRACKING 


fork=l:N-l 

Pp=phi*P*phi'+gam*Q*gam';        %  Simple  Kalman  Filter 

xp=phi*xe(:  Jc);  %  only  info,  we  have 

S=H*Pp*JT+R; 

K(:,k)=Pp*H'*inv(S); 

P=(I-K(:,k)*H)*Pp; 

PP(l,k)=P(l,l); 

PP(2,k)=P(l,2); 

PP(3,k)=P(2,l); 
PP(4,k)=P(2,2); 
inov(k+ 1  )=z(k+ 1  )-H*xp; 
xe(:,k+l)=xp+K(:,k)*inov(k+l); 


Ppm=phi*Pm*phi'+gam*Q*gam';      %  Kalman  Filter  with  Input 

xpm=phi*xem(:,k)+gam*u(k);       %  can't  use  for  real  world 

Sm=H*Ppm*H*+R; 

Km(:,k)=Ppm*H**inv(Sm); 

Pm=(I-Km(:,k)*H)*Ppm; 

PPm(l,k)=Pm(l,l); 

PPm(2,k)=Pm(l,2); 

PPm(3,k)=Pm(2,l); 

PPm(4,k)=Pm(2,2); 

inovm(k+l)=z(k+l)-H*xpm; 

xem(:,k+l)=xpm+Km(:,k)*inovm(k+l); 
end 
eval(['save  ',fname,'_trk  xem  xe  x;']); 


%  PLOTTING 
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j=l:N; 

plot(j,x(l,:)j,xe(l,:),j,xem(l,:)),grid,pause 
plot(j,x(2,:),j,xe(2,:),j,xem(2,:)),grid,pause 
plot(PP(l,:)),grid,title('Covariance  of  Pll'),pause 
plot(PP(2,:)),grid,title('Covariance  of  P12'),pause 
plot(PP(3,:)),grid,title('Covariance  of  P21'),pause 
plot(PP(4,:)),grid,title('Covariance  of  P22'),pause 
plot(K(l,:)),grid,title(*Filter  Gain  of  Kll'),pause 
plot(K(2,:)),grid,title(*Filter  Gain  of  K12'),pause 
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3.     ELEVDETEC1.M 

%  File  Name  :  elevdetecl.m 

% 

%  Purposes  : 

%     1 .  Apply  input  estimation  algorithm 

%    2.  Compute  the  input  estimates 

%    3.  Elevation  angle  only. 

% 

%  By  :  Meng,  Hsing-Han 

clear 

clc,disp('         '); 

fname=input('Give  the  name  of  target  to  be  processed  ?  \'s'); 

tgmodl=input('Maneuver  or  unmaneuver  target  ?  [m]/[u]  ','s'); 

s=input('How  many  points  of  measurement  want  to  use  ?  '); 

if  tgmodl=='m' 

eval(['load  ,,fname,'_m']); 

x=x_manu; 

z=z_manu; 
else 

eval(['load  ',fname,'_u']); 
end 

%  GET  RID  OFF  THE  FIRST  FEW  TRANSIENT  POINTS 
cpt=10; 

%  TARGET  TRACKING 

var_we=0. 1 ; 

var_e=2;  var_ed=3; 

T=l; 

phi=[l  T;0  1]; 

H=[l  0]; 
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gam=[(T*2)/(2*  1 000);T/(  1 000)] ; 

I=eye(2); 

Q=var_we; 

R=var_e; 

xe(:,l)=[z(l);0]; 

P=diag([var_e  var_ed]); 

i=0;j=0; 

fork=l:N-l 

Pp=phi*P*phi'+gam*Q*gam';        %  Simple  Kalman  Filter 

xp=phi*xe(:Jc);  %  only  info,  we  have 

w(k)=H*Pp*H'+R; 

K(:,k)=Pp*H**inv(w(k)); 

P=(I-K(:,k)*H)*Pp; 

inov(k+l  )=z(k+ 1  )-H*xp; 

xe(:,k+l)=xp+K(:,k)*inov(k+l); 

%    COMPUTE  THE  ESTIMATED  INPUT 
if  k>(cpt+s-l) 
M=0; 
for  i=l:s 

A(i,  1  )=H*(phi*M+I)*gam ; 
M=(I-K(:,k-s+i)*H)*(phi*M+I); 
Sinv(i,i)=inv(w(k-s+i)); 
end 

C=A'*Sinv; 

uest(k)=mv(C*A)*C*inov(k-s+l:k)'; 
end 
end 

%  PLOTTING 

ss=num2str(s); 

plot(uest),grid,title(['Tgtl.  ',ss,'100  measurements  to  compute  U  estimate') 
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xlabelCk'XylabeK'U  estimate'),pause 

max(abs(uest)) 

mean(uest) 

std(uest) 

%  SAVE  TO  WORKSPACE 

if  tgmodl=='m' 

eval(['save  ',fname,ss,'_est']); 
else 

eval(['save  '^namcss/u.est']); 
end 
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4.     ELEVDETEC.M 

%  File  Name  :  elevdetec.m 

% 

%  Purposes  : 

%     1.  Apply  input  estimation  algorithm 

%    2.  Compute  the  input  estimates 

%     3.  Perform  the  detection  algorithm 

%    4.  Elevation  angle  only. 

% 

%  By  :  Meng,  Hsing-Han 

clear 

clc,disp('        '); 

fname=input('Give  the  name  of  target  to  be  processed  ?  ','s'); 

s=input('How  many  points  of  measurement  want  to  use  ?  '); 

eval(['load  ,,fname,'_m']); 

x=x_manu; 

z=z_manu; 

%  GIVE  THE  THRESHOLD  VALUE 
thold=3.89; 

%  GET  RID  OFF  THE  FIRST  FEW  TRANSIENT  POINTS 
cpt=10; 

%  TARGET  TRACKING 

var_we=0. 1 ; 

var_e=2;  var_ed=3; 

T=l; 

phi=[lT;01]; 

H=[l  0]; 

gam=[(TA2)/(2*  1 000);T/(  1 000)] ; 
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I=eye(2); 

Q=var_we; 

R=var_e; 

xe(:,l)=[z(l);0]; 

P=diag([var_e  var_ed]); 

i=0;j=0; 


fork=l:N-l 

Pp=phi*P*phi'+gam*Q*gam';        %  Normal  Kalman  Filter 

xp=phi*xe(:  ,k);  %  only  info,  we  have 

w(k)=H*Pp*H'+R; 

K(:Jc)=Pp*H'*inv(w(k)); 

P=(I-K(:,k)*H)*Pp; 

PP(l,k)=P(l,l); 

PP(2,k)=P(l,2); 

PP(3,k)=P(2,l); 
PP(4,k)=P(2,2); 

inov(k+l  )=z(k+ 1  )-H*xp; 
xe(:,k+l  )=xp+K(:,k)*inov(k+l ); 


%    COMPUTE  THE  ESTIMATED  INPUT 
if  k>(cpt+s-l) 
M=0; 
fori=l:s 

A(i,l)=H*(phi*M+I)*gam; 

M=(I-K(:,k-s+i)*H)*(phi*M+I); 

Sinv(i,i)=inv(w(k-s+i)); 
end 

C=A'*Sinv; 
uest(k)=inv(C*A)*C*inov(k-s+l:k)'; 

%      DETECTION  ALGORITHM 
if  abs(uest(k))<=thold 
detec(k)=0; 
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else 

detec(k)=3; 
end 
end 
end 

m  vprd=f ind(deteoO) ; 
strpt=mvprd(l) 

%  PLOTTING 

plot(detec),grid,pause 
plot(uest(l:330)),grid 
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5.     IE_SENSITIVE.M 

%  File  Name  :  EE_sensitive.m 

% 

%  Purposes  : 

%     1.  Apply  input  estimation  algorithm. 

%    2.  Compute  the  input  estimates  from  false  maneuver  start  time 

%       using  variable  length  window. 

%    3.  Elevation  angle  only. 

%     4.  Run  this  program  for  20  times. 

% 

%  By  :  Meng,  Hsing-Han 

diary  IE3_dat315 
fort=l:20 
clear 

%    FALSE  MANEUVER  START  TIME 
MPF=315; 

%    POINTS  OF  MEASUREMENT  USED  TO  COMPUTE  INPUT  ESTIMATES 

s=50; 

%    INITIAL  CONDITION 
N=600; 
MP=300; 
amp=3; 
var_we=0.1; 
var_e=2;  var_ed=3; 
rand('normar); 
n=l:N; 

w=l*sqrt(var_we)*rand(n); 
v=sqrt(var_e)*rand(n); 
u=[zeros(  1  ,MP- 1 )  amp*ones(  1  ,N-(MP- 1 ))] ; 
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%    GET  RID  OFF  THE  FIRST  FEW  TRANSIENT  POINTS 
cptl=10; 

%    GENERATE  MANEUVERING  TARGET 
phi=[lT;01]; 
H=[l  0]; 
gam=[(TA2)/(2*  1 000);T/(  1 000)] ; 

x(:,l)=[23;0]; 
z(l)=H*x(:,l)+v(l); 

fork=l:N-l 

x(:,k+l)=phi*x(:,k)+gam*(w(k)+u(k)); 
z(k+l)=H*x(:,k+l)+v(k+l); 
end 

%    TARGET  TRACKING 
T=l; 

xe(:,l)=[z(l);0]; 
I=eye(2); 
Q=var_we; 
R=var_e; 

P=diag([var_e  var_ed]); 
i=0;j=0;M=0; 

fork=l:N-l 

Pp=phi*P*phi'+gam*Q*gam';       %  Normal  Kalman  Filter 

xp=phi*xe(:,k);  %  only  info,  we  have 

w(k)=H*Pp*IT+R; 

K(:,k)=Pp*H'*inv(w(k)); 

P=(I-K(:,k)*H)*Pp; 

PP(l,k)=P(l,l); 

PP(2,k)=P(l,2); 

PP(3,k)=P(2,l); 
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PP(4,k)=P(2,2); 
inov(k+l)=z(k+l)-H*xp; 
xe(:,k+l)=xp+K(:,k)*inov(k+l); 
ifk>MPF 
i=i+l; 

A(i,  1  )=H*(phi*M+I)*gam ; 
M=(I-K(:,k)*H)*(phi*M+I); 
Sinv(i,i)=inv(w(k)); 
C=A'*Sinv; 

uest(k)=inv(C*  A)*C*inov(MPF+ 1  :k)'; 
end 
end 

umean=mean(uest(399:N- 1 )) 
end 
diary  off 
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6.    ELEVESDCR.M 

%  File  Name  :  elevesdcr.m 

% 

%  Purposes  : 

%     1.  Apply  input  estimation  algorithm 

%     2.  Compute  the  input  estimates 

%     3.  Perform  the  detection  algorithm 

%     4.  Apply  the  correction  algorithm  when  the  absolute  value  of 

%       input  estimates  less  than  or  equal  to  the  threshold  value 

%     5.  Elevation  angle  only. 

% 

%  By  :  Meng,  Hsing-Han 

clear 

clc,disp('        '); 

fname=input('Give  the  name  of  target  to  be  processed  ?  ','s'); 

s=input('How  many  points  of  measurement  want  to  use  ?  '); 

eval(['load  ',fname,'_m']); 

x=x_manu; 

z=z_manu; 

%  GIVE  THE  THRESHOLD  VALUE 
th=3.89; 

%  GET  RID  OFF  THE  FIRST  FEW  TRANSIENT  POINTS 
cpt=10; 

%  TARGET  TRACKING 

var_we=0. 1 ; 
var_e=2;  var_ed=3; 
T=l; 
phi=[lT;0  1]; 
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H=[l  0]; 

gam=[(TA2)/(2*  1 000);T/(  1 000)] ; 

I=eye(2); 

Q=var_we; 

R=var_e; 

xe(:,l)=[z(l);0]; 

P=diag([var_e  var_ed]); 

i=0;j=0; 


fork=l:N-l 

Pp=phi*P*phi'+gam*Q*gam';        %  Normal  Kalman  Filter 

xp=phi*xe(:,k);  %  only  info,  we  have 

w(k)=H*Pp*H'+R; 

K(:*)=Pp*H**inv(w(k)); 

P=(I-K(:,k)*H)*Pp; 

PP(l,k)=P(l,l); 

PP(2,k)=P(l,2); 

PP(3,k)=P(2,l); 

PP(4,k)=P(2,2); 

inov(k+ 1  )=z(k+ 1  )-H*xp; 

xe(:,k+l)=xp+K(:,k)*inov(k+l); 


%    COMPUTE  THE  ESTIMATED  INPUT 
if  k>(cpt+s-l) 
M=0; 
for  i=l:s 

A(i,l)=H*(phi*M+I)*gam; 

M=(I-K(:,k-s+i)*H)*(phi*M+I); 

Sinv(i,i)=inv(w(k-s+i)); 
end 

C=A'*Sinv; 
uest(k-s+l)=inv(C*A)*C*inov(k-s+l:k)'; 

%      DETECTION  AND  CORRECTION  ALGORITHM 
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if  abs(uest(k-s+l))<=th 
xe(:,k+l)=xe(:,k+l)+M*gam*uest(k-s+l); 
P=P+(M*gam)*inv(C*A)*(M*gam)'; 
end 
end 
end 
eval(['save  ',fname,'  uest  xe  inov  PP  K;']); 

%  PLOTTING 

j=l:600; 

plot(uest),grid,title('Estimated  Input'),pause 

plot(j,x(l,:)J,xe(l,:)),grid 

title('Real  Elevation  Angle  vs.  Estimated  Elevation  Angle'),pause 

plot(j,x(2,:),j,xe(2,:)),grid 

title('Real    Elevation    Angle    Derivate    vs.    Estimated    Elevation   Angle 

Derivate'),pause 

plot(PP(l  ,:)),grid,pause 

plot(PP(2,:)),grid,pause 

plot(PP(3,:)),grid,pause 

plot(PP(4,:)),grid,pause 

plot(K(l  ,:)),grid,pause 

plot(K(2,:)),grid 
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7.     ELEVESDCRA.M 

%  File  Name  :  elevesdcr_a.m 

% 

%  Purposes  : 

%     1.  Apply  input  estimation  algorithm 

%     2.  Compute  the  input  estimates 

%     3.  Perform  the  detection  algorithm 

%    4.  Apply  the  correction  algorithm  when  the  absolute  value  of 

%       input  estimates  greater  than  the  threshold  value 

%     5.  Elevation  angle  only. 

% 

%  By  :  Meng,  Hsing-Han 

clear 

clc,disp('        '); 

fname=input('Give  the  name  of  target  to  be  processed  ?  ','s'); 

s=input('How  many  points  of  measurement  want  to  use  ?  '); 

eval(['load  \fname,'_m']); 

x=x_manu; 

z=z_manu; 

%  GIVE  THE  THRESHOLD  VALUE 
thold=3.89; 

%  GET  RID  OFF  THE  FIRST  FEW  TRANSIENT  POINTS 
cpt=10; 

%  TARGET  TRACKING 

var_we=0. 1 ; 
var_e=2;  var_ed=3; 
T=l; 
phi=[lT;01]; 
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H=[l  0]; 

gam=[(TA2)/(2*1000);T/(1000)]; 

I=eye(2); 

Q=var_we; 

R=var_e; 

xe(:,l)=[z(l);0]; 

P=diag([var_e  var_ed]); 

i=0;j=0; 


fork=l:N-l 

Pp=phi*P*phi'+gam*Q*gam';        %  Normal  Kalman  Filter 

xp=phi*xe(:,k);  %  only  info,  we  have 

w(k)=H*Pp*H'+R; 

K(:,k)=Pp*H'*inv(w(k)); 

P=(I-K(:,k)*H)*Pp; 

PP(l,k)=P(l,l); 
PP(2,k)=P(l,2); 

PP(3,k)=P(2,l); 
PP(4,k)=P(2,2); 
inov(k+ 1  )=z(k+ 1  )-H*xp; 
xe(:,k+l)=xp+K(:,k)*inov(k+l); 


%    COMPUTE  THE  ESTIMATED  INPUT 
ifk>(cpt+s-l) 
M=0; 
fori=l:s 

A(i,l)=H*(phi*M+I)*gam; 

M=(I-K(:,k-s+i)*H)*(phi*M+I); 

Sinv(i,i)=inv(w(k-s+i)); 
end 

C=A'*Sinv; 
uest(k-s+l)=inv(C*A)*C*inov(k-s+l:k)'; 

%      DETECTION  AND  CORRECTION  ALGORITHM 
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if  abs(uest(k-s+l))>=thold 
xe(:,k+l)=xe(:,k+l)+M*gam*uest(k-s+l); 
P=P+(M*gam)*inv(C*A)*(M*gam)'; 
end 
end 
end 
eval(['save  \fname,'a  uest  xe  inov  PP  K;']); 

%  PLOTTING 

j=l:N; 

plot(uest),grid,title('Estimated  Input'),pause 
plot(j,x(l,:),j,xe(l,:)),grid 

title('Real  Elevation  Angle  vs.  Estimated  Elevation  Angle'),pause 

plot(j,x(2,:),j,xe(2,:)),grid 

title('Real    Elevation    Angle    Derivate    vs.    Estimated    Elevation   Angle 

Derivate'),pause 

plot(PP(  1  ,:)),grid,pause 

plot(PP(2,:)),grid,pause 

plot(PP(3,:)),grid,pause 

plot(PP(4,:)),grid,pause 

plot(K(  1  ,:)),grid,pause 

plot(K(2,:)),grid 
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8.     ELEVESDCR_D.M 

%  File  Name  :  elevesdcr_d.m 

% 

%  Purposes : 

%     1 .  Apply  input  estimation  algorithm 

%     2.  Compute  the  input  estimates  using  fixed  length  window 

%     3.  Use  the  computed  input  estimations  to  elevesdcr_d  1  .m 

%     4.  Elevation  angle  only. 

% 

%  By  :  Meng,  Hsing-Han 

clear 

clc,disp('        '); 

fname=input('Give  the  name  of  target  to  be  processed  ?  ','s'); 

s=input('How  many  points  of  measurement  want  to  use  ?  '); 

eval(['load  ',fname,'_m']); 

x=x_manu; 

z=z_manu; 

%  GET  RID  OFF  THE  FIRST  FEW  TRANSIENT  POINTS 
cpt=10; 

%  TARGET  TRACKING 

var_we=0. 1 ; 

var_e=2;  var_ed=3; 

T=l; 

phi=[lT;01]; 

H=[l  0]; 

gam=[(T^2)/(2*  1 000);T/(  1 000)] ; 

I=eye(2); 

Q=var_we; 

R=var_e; 
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xe(:,l)=[z(l);0]; 
P=diag([var_e  var_ed]); 
i=0;j=0; 


fork=l:N-l 

Pp=phi*P*phi'+gam*Q*gam';        %  Normal  Kalman  Filter 

xp=phi*xe(:,k);  %  only  info,  we  have 

w(k)=H*Pp*H*+R; 

K(:Jk)=Pp*H*inv(w(k)); 

P=(I-K(:,k)*H)*Pp; 

PP(l,k)=P(l,l); 

PP(2,k)=P(l,2); 

PP(3,k)=P(2,l); 

PP(4,k)=P(2,2); 

inov(k+ 1  )=z(k+ 1  )-H*xp; 

xe(:,k+l)=xp+K(:,k)*inov(k+l); 


%    COMPUTE  THE  ESTIMATED  INPUT 
if  k>(cpt+s-l) 
M=0; 
fori=l:s 

A(i,l)=H*(phi*M+I)*gam; 
M=(I-K(:,k-s+i)*H)*(phi*M+I); 
Sinv(i,i)=inv(w(k-s+i)); 
end 

C=A'*Sinv; 

uest(k)=inv(C*A)*C*inov(k-s+l:k)'; 
end 
end 

ss=num2str(s); 

eval(['save  '.fname/J.ss.'uest  uest;']); 

elevesdcr  dl 
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9.     ELEVESDCR_D1.M 

%  File  Name  :  elevesdcr_dl.m 

% 

%  Purposes  : 

%     1 .  Use  the  computed  input  estimates  to  the  Kalman  filter 

%    2.  Elevation  angle  only. 

% 

%  By  :  Meng,  Hsing-Han 

clear 

clc,disp('        '); 

fname=input('Give  the  name  of  target  to  be  processed  ?  ','s'); 

s=input('How  many  points  of  measurement  want  to  use  ?  '); 

amp=input('Give  the  magnitude  of  input  : '); 

ss=num2str(s); 

eval(['load  ,,fname,,_,,ss,,uest;']); 

eval(['load  ',fname,'_m']); 

x=x_manu; 

z=z_manu; 

dd=amp/mean(uest(600:999)) 

%  TARGET  TRACKING 

var_we=0.1; 

var_e=2;  var_ed=3; 

T=l; 

phi=[lT;01]; 

H=[l  0]; 

gam=[(TA2)/(2*  1 000);T/(  1 000)] ; 

I=eye(2); 

Q=var_we; 

R=var_e; 

xe(:,l)=[z(l);0]; 
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P=diag([var_e  var_ed]); 
i=0;j=0; 

fork=l:N-l 

Pp=phi*P*phi'+gam*Q*gam'; 

xp=phi*xe(:,k)+gam*uest(k); 

w(k)=H*Pp*H'+R; 

K(:,k)=Pp*H,*inv(w(k)); 

P=(I-K(:,k)*H)*Pp; 

PP(l,k)=P(l,l); 

PP(2,k)=P(l,2); 

PP(3,k)=P(2,l); 

PP(4,k)=P(2,2); 

inov(k+ 1  )=z(k+ 1  )-H*xp; 

xe(:,k+l)=xp+K(:,k)*inov(k+l); 
end 

eval(['save  ',fname,ss,'d4  uest  xe  inov  PP  K;']); 

%  PLOTTING 

j=l:N; 

plot(uest),grid,title('Estimated  Input'),pause 

plot(j,x(l,:),j,xe(l,:)),grid 

title('Real  Elevation  Angle  vs.  Estimated  Elevation  Angle'),pause 

plot(j,x(2,:),j,xe(2,:)),grid 

title('Real    Elevation    Angle    Derivate    vs.    Estimated    Elevation   Angle 

Derivate'),pause 

plot(PP(  1  ,:)),grid,pause 

plot(PP(2,:)),grid,pause 

plot(PP(3,:)),grid,pause 

plot(PP(4,:)),grid,pause 

plot(K(  1  ,:)),grid,pause 

plot(K(2,:)),grid 
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